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I/"") ' Abstract. In these proceedings we present recent efforts to understand the energetics of magneto- 

hydrodynamic (MHD) turbulence driven by the magneto-rotational instability (MRI). These studies 
are carried out in the local (shearing box) approximation using the Athena simulation code. Athena 
is a higher order Godunov algorithm based on the piecewise parabolic method (PPM), the cor- 
ner transport upwind (CTU) integration algorithm, and the constrained transport (CT) algorithm 
for evolving the magnetic field. This algorithm is particularly suited for these studies owing to the 
conservation properties of a Godunov scheme and the particular implementation of the shearing box 
^-j- . source terms used here. We present a variety of calculations which may be compared directly to pre- 

^vq ' viously published results and discuss them in some detail. The only significant discrepancy found 

between the results presented here and in the published literature involves the turbulent heating rate. 
We observe the presence of recurrent channel solutions in calculations involving a mean vertical 
magnetic field and the associated time lag between the energy injection and thermalization rate. We 
ON , also present the results of a shearing box calculation which includes an optically thin radiative term 

with a cooling rate selected to match the turbulent heating rate. Some properties which we find uni- 
formly present in all of the calculations presented here are compressible fluctuations, spiral waves 
and weak shocks. It is found that these compressible modes dominate the temporal fluctuations in 
the probability distribution functions for most of the thermodynamic variables; only the specific 
| entropy is relatively immune to their effects. 
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INTRODUCTION 

The last 10 years have borne witness to a substantial improvement in our understanding 
of the role of sub-thermal magnetic fields in accretion disk physics [1]. The magneto- 
rotational instability (MRI) plays the central role in driving magnetohydrodynamic 
(MHD) turbulence which in turn transports mass and angular momentum. A large body 
of work now exists in the literature which has elucidated a great many aspects of the 
MRI. Yet surprisingly, one aspect of accretion disk physics has yet to receive significant 
attention; this is the role of energetics 1000]. 

One of the simplest energetics questions which could be addressed with numerical 
simulations is, what energy reservoir is tapped, how is the energy distributed amongst 
its kinetic and magnetic forms, and where does it ultimately go. Some elements of this 
question have been answered repeatedly by a number of researchers. The differential 
rotation is the ultimate energy source which is tapped via the MRI and injected into 
the kinetic and magnetic energy. The energy stored in the form of magnetic and kinetic 
energy may be exchanged yD, but ultimately, it is dissipated and ends up as thermal 
heating. Interestingly enough, this last stage is relatively unimportant as far as the 
MRI itself is concerned since there is only a weak pressure dependence on the growth 
rates |0] and saturation amplitudes |0]. Nevertheless, it is the last stage of this process 



which results in radiation and allows a comparison between theoretical calculations and 
observations. 

There are two circumstances under which the thermal state of the plasma can be easily 
understood to have a significant impact on an MRI turbulent accretion disk. The first 
is in determining the vertical structure of a stratified accretion disk. The scale height 
of an accretion disk is determined by balancing the thermal pressure gradients against 
gravity. The disk scale height is thus a measure of the mean or central gas pressure, and 
is determined directly from the balance of turbulent heating (as well as other potential 
sources) and radiative losses. (This is a question we intend to explore in the near future.) 
The second is by the strong temperature dependence on the ionization state, resistivity, 
and other kinetic effects which directly modify the induction equation such as ambipolar 
diffusion or the Hall conductivity. For example, the degree to which the MRI plays a 
role in protostellar accretion will depend heavily on these factors It has also been 

suggested that the temperature sensitivity of the resistivity could serve as a switch which 
regulates the outburst activity in dwarf nova disks II 1 OL I X 111 . 

The focus of this paper is twofold with the ultimate goal to seek a better understanding 
of the energetics in MRI turbulent accretion disks. This goal is ultimately achievable as a 
result of the recent development of a fully second order accurate, conservative Godunov 
algorithm for ideal MHD [12|]. Using the shearing box formalism, one is assured that any 
energy injected into the computational domain ultimately results in thermally heating 
the plasma without directly modeling the physical dissipation processes. Nevertheless, 
since this is a new computational algorithm, and is the first time that a fully second order 
Godunov algorithm has been applied to the study of the MRI we are compelled to repeat 
a series of calculations which can be directly compared with published results. The 
second goal of this paper is to validate this new computational algorithm when applied 
to the study of MRI driven MHD turbulence in the shearing box and to understand 
any differences which arise between the results presented here and published in the 
literature. To that end, we present results for adiabatic calculations with three different 
field geometries and a calculation including radiative cooling. 

NUMERICAL CALCULATIONS 
Shearing Box Approximation 

The studies presented here make use of the local "shearing box" formalism ifoll . In 
this approach, one focuses on a local section of an accretion disk at some fiducial radius 
Rq and in a reference frame co-rotating at the local angular velocity Q.q. Expanding 
the ideal MHD system of equations in cylindrical coordinates for small deviations 
8r/Ro <C 1, etc. and relabeling the coordinates (dr,R$d<l) ,dz) as {dx,dy,dz) one obtains 
the ideal MHD system in a Cartesian domain plus and additional set of source terms 
which represent the Coriolis and tidal gravity forces. For a Keplerian accretion disk one 
obtains 



3T + V.(pv) = 



(1) 



pQ.l(3xi-zk) -2Q k x pv (2) 

(3) 

Q. 2 p\-(3xi-zk) (4) 

where p is the mass density, pv is the momentum density, B is the magnetic field, and 
E is the total energy density. The total pressure P* = P + (B • B)/2 where P is the gas 
pressure, and the total energy density E is related to the internal energy density e via 

£ = £ + p(vv)/2 + (B-B)/2. (5) 

This system of equations is closed with an equation of state relating the pressure to the 
density and internal energy. We use an ideal gas equation of state for which P = (7— l)e, 
and take 7 = 5/3. In the calculations presented in this paper we further simplify this 
system of equations by neglecting the vertical z-component of gravity. As a result, the 
calculations presented here are applicable to the central ~ 1 scale height of the accretion 
disk. 

We present calculations both with and without radiative cooling. In those calculations 
which include optically thin radiative cooling we choose a simple form for the cooling 
term which is appropriate for bremsstrahlung radiation. In these cases we add a term to 
the right hand side of equation © of the form 

-Ap 2 (e/p) 1 / 2 (6) 
and choose the coefficient A such that the cooling rate matches the turbulent heating rate. 

Numerical Method 

The numerical method we use to solve the MHD equations in the shearing box ap- 
proximation can be described as having three basic parts. These are the basic integration 
algorithm, the treatment of source terms within this algorithm and the application of the 
boundary conditions. We briefly describe each of these in turn in what follows. 

Integration Algorithm 

We solve the shearing box system of MHD equations using the recently developed 
simulation code, Athena. This code was constructed to solve the system of ideal hy- 
drodynamics (HD) and magnetohydrodynamics (MHD) by combining a few key algo- 
rithms. The integration algorithm is based upon the directionally unsplit corner transport 
upwind (CTU) integration algorithm of Collela lfl4l[l5ll . The CTU integration algorithm 
can be formally understood as a procedure for correcting the interface states in the piece- 
wise parabolic method (PPM) HI 611 to account for multidimensional wave propagation. 



dpv 



+ V-(pvv-BB) + VP* 



dB 



+ V-(vB-Bv) 



dE 
~dt 



+ V-((E + P*)v-B(B-v)) 



We have combined this algorithm with the method of constrained transport (CT) ifnll 
to evolve the magnetic field and explicitly preserve the V ■ B = condition. A detailed 
description of the construction of this integration algorithm for MHD in two dimensions 
is in press lfl2ll and a description of the three dimensional algorithm is in preparation. 

At the time of this writing, detailed code verification calculations, user and program- 
mer's guides, download information, etc. can be found on the Athena home page at 
|http : / / www .astro . prince ton . edu/~ j stone/ athena/ |, See also the arti- 
cle by J. Stone in these proceedings for more information on the Athena integration 
algorithm and a description of the first application results. 



Source Terms 

One of our principal goals in this work is to study the energetics in MRI turbulent 
accretion disks. Hence, the evolution of the energy equation © is of particular interest. 
Energy conservation can be restored to this system of equations by solving for the evo- 
lution of the total energy E t = (E + p4>), including the tidal potential, 4> = — 1.5£lfyc 2 . 
With this choice, one is guaranteed by the conservative nature of the equations, that en- 
ergy can be exchanged between its kinetic, magnetic, and thermal forms, but cannot be 
lost in the form of truncation error. The only route by which energy can be added to or 
removed from the computational domain is through the boundaries. In fact, making use 
of the shearing periodic boundary conditions one can show lfl3ll that the volume average 
of the total energy < E + p$ > evolves as 

|- < E + pcl> >= f (pv x 8v y - B x B y ) dydz (7) 

at L y L z Jx 

where 8v y = v y + 1.5Q.QX is the angular velocity fluctuation, and the integral is taken 
over one of the bounding x-faces. This is precisely the route by which the MRI taps the 
kinetic energy in the form of differential rotation in the shearing box and amplifies the 
magnetic field. In the calculations presented here, the integral relation in equation Q is 
satisfied to numerical roundoff error. 

We have found that the treatment of the source terms in the momentum equation © 
is particularly important for shearing box calculations with Athena. Note that unlike 
the energy equation, the momentum equation cannot be recast in a conservative form. 
Moreover, there is a simple limit in which the x- and y-momentum equations are strongly 
coupled via the source terms and therefore should not be evolved independently. This 
strong coupling limit describes the solution for uniform epicyclic oscillation modes. For 
epicyclic oscillations there is a conserved kinetic energy 

E epi = l -p{v 2 x + A{8v y ) 2 ) . (8) 

We have found that implementing the source terms in Athena in such a way that the 
epicyclic kinetic energy is conserved is extremely important. 

In order to clarify why the conservation of the epicyclic kinetic energy should have 
such a large impact on the solution, note that it is the long wavelength modes which 



contain a dominant share of the energy. This is a natural consequence of the shearing 
box formalism in which the driving scale is effectively equal to the radial size of the 
computational domain and the dissipation which occurs in the form of a turbulent cas- 
cade. In addition to long wavelength modes which fit into the computational domain, the 
shearing box MHD equations (and boundary conditions) also support uniform epicyclic 
oscillation modes. It is reasonable to expect that some coupling between long wave- 
length modes (which fit in the computational domain) and epicyclic oscillation modes 
may result from numerical truncation error. By conserving the epicyclic kinetic energy, 
one ensures that this coupling can not set up a positive feedback loop amplifying the 
kinetic energy in the box. 

In order to conserve the epicyclic kinetic energy we evolve not the v-momentum equa- 
tion, but the equivalent conservation law for the angular momentum fluctuation p8v y . 
When written in this form, uniform epicyclic motion is consistent with the vanishing of 
the flux gradient terms. This leaves the x- and y-momentum fluctuations evolution to be 
controlled by the resulting source terms. One can show that evaluating the source terms 
using Crank-Nicholson preserves the epicyclic kinetic energy. 

For those calculations which include radiative cooling, we apply the operator splitting 
technique and split the radiative cooling equation 

§ = -Ap 2 ( £ /p) 1 /2 (9 ) 

from the rest of the shearing box evolution equations and solve for the radiative cooling 
at constant density. 



Boundary Conditions 

The boundary conditions which we employ in these calculations are strictly periodic 
in the y- and z-directions and shearing periodic in the ^-direction. One result of these 
boundary conditions is that so long as the volume averaged ^-component of the magnetic 
field is zero, the mean flux through the computational domain is fixed in time. The 
implementation of the boundary conditions in this work preserves this condition to 
roundoff error for the x- and z-components of the magnetic field and to truncation error 
for the y-component of the magnetic field. Moreover, the remapping of the magnetic 
field in the ghost zones at the ^-boundaries of the computational domain preserves 
the V • B = condition to roundoff error. The details of how this is achieved will be 
described in a future paper. 

RESULTS 

In this section we will present results from three local shearing box simulations with 
different magnetic field geometries. The cases which we consider here include a mean 
toroidal field, a mean vertical field, and a zero net- flux case. Our choice of computational 
domain and initial conditions mirrors that used in [1 1 3fl . Namely, we choose a computa- 
tional domain given by —0.5 < x < 0.5, —it < y < K, —0.5 < z < 0.5 and resolve it on 
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FIGURE 1. Time evolution of the volume averaged magnetic, kinetic, total and thermal energy for the 
mean toroidal field calculation Y2. Time is in units of the orbital period. 




FIGURE 2. Time evolution of the volume averaged Maxwell and Reynolds stress for the mean toroidal 
field calculation Y2. Time is in units of the orbital period. 



a 64 x 128 x 64 grid, except where stated otherwise. The initial velocity components 
v x = v z = and v y = -1.5Q.qx. We also choose = 10~ 3 , the initial density p = 1 
and pressure P = 10 > 6 . To these initial conditions we add 0.1% white noise adiabatic 
density and pressure perturbations. For reference, in the following analysis we choose 
the zero point of the tidal potential <J» = (0.125 — 1.5x 2 )Q.q so that initially the mean 
tidal potential energy < p4> >= 0. 



Uniform Toroidal Magnetic Field 

The first case we consider has an initial magnetic field given by B x = B z = and 
By = const, with /3 = 2P/B 2 = 100. This model calculation can be compared directly to 
models Yl and Yll in [1 1 3fl . In figure we present the time history of the volume av- 
eraged Maxwell and Reynolds stress. The volume averaged Maxwell stress < —B x B y > 
/Pq and Reynolds stress < pv x 8v y > /Pq have a mean value of 0.040 ± 8.0 x 10~ 3 and 
8.7 x 10~ 3 ±3.4 x 10 3 respectively where the error is given by one standard deviation. 
The ratio of the Maxwell to Reynolds stress is ~ 4.4. These values are in good agree- 
ment with the results of 11311 . Similarly we find that the saturation value for the kinetic 
and magnetic energy density shown in figure dare in good agreement with the results of 

Q. _ 

The results presented here differ from those in H13I1 when we consider the evolution of 
the thermal energy density. In figure [l] we present the evolution of the volume averaged 
thermal energy density < £ > and the total energy density < E + p4> >. As noted 




FIGURE 3. Time evolution of the volume averaged magnetic, kinetic, total and thermal energy for the 
mean vertical field calculation Z2. Time is in units of the orbital period. 
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FIGURE 4. Time evolution of the volume averaged Maxwell and Reynolds stress for the mean vertical 
field calculation Z2. Time is in units of the orbital period. 



previously, the total energy in the shearing box is conserved and its volume average 
can increase only through the differences of the flux across the radial boundaries. From 
figure [l] we see that the internal and total energies track one another very closely as 
they must after the kinetic and magnetic energy densities reach saturation. Applying a 
linear curve fit to the evolution of the total energy density after saturation and applying 
equation (Q, one may extract the surface averaged value of the Maxwell plus Reynolds 
stress one finds a very good agreement with the volume averages noted above. When 
compared to the results in H 1 3fl we find that the turbulent heating rate is increased by 
nearly an order of magnitude. This difference is likely attributable to the lack of a 
resistive heating term in the evolutionary equation for the internal energy. 



Uniform Vertical Magnetic Field 

Next we consider the case where the magnetic field is initialized with B x = B y = and 
B z = const, with /3 = 2P/B 2 = 400. This model calculation can be compared directly 
to models Z4, Z7, Z19 and Z22 in 11311 . In figure |3] we present the temporal evolution 
of the volume averaged kinetic, magnetic, thermal and total energies, and in figure|Hthe 
Maxwell and Reynolds stress. It is immediately apparent upon inspection that this field 
configuration is much more time variable as a result of recurrent channel flows [18]. The 
volume averaged Maxwell stress < —B x B y > /Pq and Reynolds stress < pv x 8v y > /Pq 
have a mean value of 0.33 ± 0. 1 8 and 0.05 1 ± 0.040 respectively where the error is given 
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FIGURE 5. Time evolution of the time derivative of the volume averaged total energy and thermal en- 
ergy. This figure demonstrates a clearly discernible delay between the energy injection and thermalization 
rate. Time is in units of the orbital period. 



by one standard deviation. These values are approximately 4 times larger than for the 
mean toroidal field configuration. The ratio of the Maxwell to Reynolds stress is 6.4. 
These numbers are in good agreement with II 1 311 . 

The evolution of the volume averaged internal energy < £ > in figure |3] shows 
approximately a factor of 7 times larger heating rate than was found in JH. It is 
interesting to note that the time variability of the Maxwell and Reynolds stress in these 
flows is sufficiently large that it results in a clearly discernible stair stepping on the 
evolution of the total and internal energy density. 

It was shown in [18] that the presence of recurrent channel flows results in a clearly 
observable time delay between the time derivative of the volume average of the total 
energy < E + p4> > and the time derivative of the volume average of the thermal energy 
density < £ >. Note that these quantities are equivalent to the energy injection rate and 
the energy thermalization rate. Differentiating the total and thermal energy plots in figure 
|3]and normalizing the derivative to the ratio of the initial thermal energy density and the 
orbital period, one obtains an equivalent plot, shown in figure to that presented in 
iflSll . Despite the presence of high frequency oscillations in the thermal energy evolution 
(which are a result of magnetosonic oscillations) there is a clearly discernible time lag 
between the peak in the energy injection rate and the average turbulent dissipation rate. 
The results presented in figure |5] also agree with 111 8ll in that the thermalization rate 
overshoots the energy injection rate at its maximum. The fact that the results presented 
here are for an ideal MHD fluid, while those in 111 8ll are for a magnetic Reynolds 
number ReM = 1 provides strong support for both calculations. The time lag between 
the energy injection rate and thermalization rate is a direct result of the channel flows 
and is essentially independent of the nature of the resistivity. Note that this delay is not 
observed for any other field configuration. 




FIGURE 6. Time evolution of the volume averaged magnetic, kinetic, total and thermal energy for the 
zero net field calculation S2. Time is in units of the orbital period. 




FIGURE 7. Time evolution of the volume averaged Maxwell and Reynolds stress for the zero net field 
calculation S2. Time is in units of the orbital period. 

Zero Net Flux 

Next we consider the case where the magnetic field is initialized with B x = B y = and 
B z = BQ%m(2nx/L x ) with /3 = 1PjB\ = 4000. It was shown in JH] that the saturation 
of the MRI is independent of the initial value of /3 . As a result, this model calculation 
can be compared to the results of models SZ1, SZ2 and SZ3 in lfl9ll . In figure |6] we 
present the temporal evolution of the volume averaged kinetic, magnetic, thermal and 
total energies, and in figure |7] the Maxwell and Reynolds stress. The volume averaged 
Maxwell stress < —B x B y > /Pq and Reynolds stress < pv x 8v y > JPq have a mean value 
of 9.5 x 10~ 3 ±2.8 x 10~ 3 and 2.9 x 10~ 3 ± 1.3 x 10~ 3 respectively where the error is 
given by one standard deviation. These values are approximately 4 times smaller than for 
the mean toroidal field configuration. The saturation amplitude of the magnetic energy 
density agrees well with the results of lfl9ll . Another feature of these flows is that they 
strongly excite compressional fast magnetosonic waves which propagate both inward 
and outward in the computational domain. These compressional waves are primarily 
responsible for the rapid time variability in the mean Reynolds stress as shown in figure 
|7J The presence of these waves was also noted in II 1911 . 



Compressibility Effects 



A rather dominant feature which we have observed in our shearing box simulations 
is the formation of spiral waves and weak shocks. Figure [U shows an example of a 




FIGURE 8. Grayscale image of the density in model S0C1 showing the presence of spiral shock waves. 
The density ranges from 0.8 (black) to 1.2 (white). 



weak shock in the density in the model S0C1 calculation described in the following 
section. Watching the time evolution of the density one finds the recurrent formation 
of weak spiral shocks which rapidly propagate across the grid, dissipate and at some 
later time are regenerated. One way to study this phenomena is to look at the probability 
distribution function (PDF) for the density. Perhaps the most illuminating aspect of the 
density PDF is the time evolution of the standard deviation as shown in figure |5| The 
first feature to notice in this figure is that all three models show a rapid increase in the 
density dispersion as the system makes the transition to a turbulent state followed by 
a gradual decrease toward some small positive value. Note that the larger the heating 
rate, the faster the standard deviation of the density decreases. These results support the 
conclusion that the decrease in the density dispersion as the system evolves is a simple 
consequence of the increasing thermal energy density driving the system toward the 
incompressible limit. Calculations initiated with a larger initial pressure also support 
this conclusion. Interspersed on this plot are regions of rapid variability, which are 
unresolved in the figure, but have a well defined characteristic period described by 
fast magnetosonic oscillations with an oscillation period T ~ L x /cf ast . These coherent 
magnetosonic oscillations appear to be the direct result of the localized dissipation 
of magnetic energy in current sheets via "numerical resistivity" and the shearing box 
boundary conditions which effectively make the shearing box a resonant cavity like 
structure at wavelengths comparable to the box size. In the Y2 and Z2 models at 
late times we find that the PDF for the density, temperature, and specific entropy are 
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FIGURE 9. Time evolution of the standard deviation of the density PDF for the three adiabatic calcu- 
lations S2, Y2, and Z2. Time is in units of the orbital period. 



well described by a Gaussian distribution. At early times, the density and temperature 
PDFs for all calculations presented here are highly time variable resulting from the 
compressible modes in the shearing box. 

It is worth pointing out that the initial pressure in the models presented here was 
selected to correspond to approximately the central 1 scale height of a vertically stratified 
accretion disk. In a time steady accretion disk calculation, losses due to radiative cooling 
will balance the turbulent heating and the compressibility will be non-negligible. These 
results suggest that the localized dissipation of magnetic energy in current sheets and the 
resulting compressional waves may play an important role determining the temperature 
fluctuations and hence the spectra of accretion disks. 

Radiative Zero Net Flux 

In this section we present the results of a zero net flux calculation including radiative 
cooling. This calculation is identical to the zero net flux calculation presented earlier 
with the exception that it is resolved on a 128 x 256 x 128 grid. In figure [10] we 
present the time evolution of the volume averaged kinetic, magnetic, thermal and total 
energies. The radiative cooling rate was selected in this calculation to match the turbulent 
heating rate in the saturated state. As a result the total and thermal energy density show 
little time evolution. Note also the highly time variable kinetic energy density which 
is a result of the magnetosonic oscillations and spiral shock waves which are present 
in this calculation. In figure [TT] we present the volume average of the Maxwell and 
Reynolds stress. The volume averaged Maxwell stress < —B x B y > /Pq and Reynolds 
stress < pv x Sv y > /Pq have a mean value of 3.2 x 10~ 3 ± 1.7 x 10~ 4 and 1.4 x 10~ 3 ± 
6.4 x 10 -4 respectively where the error is given by one standard deviation. These values 
are approximately a half of the values found in the calculations presented earlier at half 
the resolution. This strong resolution dependence is not found in either the mean toroidal 
or mean vertical field calculations. Note that these values for the Maxwell and Reynolds 
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FIGURE 10. Time evolution of the volume averaged magnetic, kinetic, total and thermal energy for the 
zero net field calculation S0C1. Time is in units of the orbital period. 
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FIGURE 11. Time evolution of the volume averaged Maxwell and Reynolds stress and the radiative 
cooling rate for the zero net field calculation S0C1. Time is in units of the orbital period and the radiative 
cooling rate is in units of (cq/t), the initial internal energy divided by the orbital period. 



stress are not a result of the presence of radiative cooling in these calculations, as we 
have also performed non-radiative calculations at this resolution and find essentially 
identical values. 

In figure [TT] we present the volume averaged radiative cooling rate in units of (Eq/x), 
the ratio of the initial internal energy density to the orbital period. Note that, as a 
consistency check, one can relate the radiative cooling rate the the mean Maxwell 
and Reynolds stress shown in figure [TT] via equation (Q under the assumption that the 
turbulent heating is balanced by the radiative cooling. As should be expected, one finds a 
very good agreement from this test. There are two main features present in this plot, the 
rapid rise as the system makes the transition to the turbulent state, and the oscillations 
resulting from the compressional magnetosonic waves and spiral shocks. While the 
high frequency oscillations are strongly determined by the boundary conditions of the 
shearing box, it is suggestive that the coupling of compressional modes to radiative 
losses may have an impact on the observed spectrum from accretion disks. 



CONCLUSIONS 

We have performed a series of shearing box calculations for three different field ge- 
ometries with the recently developed simulation code Athena. This is the first time that 
a fully second order accurate Godunov algorithm has been applied to the study of the 
magnetorotational instability in the shearing box formalism. As such, comparing the cal- 
culations presented here to previously published results provides an independent confir- 



mation. In general we have found that for non-radiative calculations the saturation am- 
plitudes of the magnetic and kinetic energy, the Maxwell and Reynolds stress are in good 
agreement with the results published in H13LI19I1 . We also found from the non-radiative 
calculations that the energy dissipation rate is underestimated in lfl3ll and presumably 
also in H 1 911 by a factor of 7 — 10. This difference is likely attributable to the lack of a 
resistive heating term in the evolutionary equation for the internal energy. 

One of the consequences of energy conservation is the accurate treatment of the dissi- 
pation of magnetic energy into thermal. This has been found to influence the solution in 
the shearing box in two ways. The first is the secular heating of the plasma as must occur 
since over a long time average the energy injection and thermalization rates must bal- 
ance. This has the tendency to push the system toward an incompressible configuration 
in which the PDF of the thermodynamic variable tends toward a Gaussian distribution. 
The second is the generation of compressible magnetosonic waves and spiral shocks. 
This is presumably a result of the fact that in ideal MHD calculations the numerical 
resistivity is nonlinear leading to the dissipation of the magnetic field predominantly in 
the form of current sheets. This localized heating process generates compressional waves 
which effectively resonate in the shearing box. These compressional modes, which are 
most important at early times, have a strong influence on the PDF for the density and 
temperature in these simulations. In radiatively cooling calculations they also have a di- 
rect influence on the radiative cooling rate. These results are suggestive that magnetic 
field dissipation in current sheets and the resulting compressional modes may have an 
observable influence on an accretion disk spectrum. 

We have also observed the recurrent channel solutions in simulations with a mean 
vertical magnetic field and the associated time lag between the energy injection rate and 
energy thermalization rate described in [18]. The fact that the calculations presented 
here were performed with ideal MHD, while those in 111 811 have a magnetic Reynolds 
number Re^ = 1 provides strong support that these results are essentially independent 
of the magnetic resistivity, so long as the MRI can operate. 

The energetics of turbulent accretion disks remains one of the least well explored and 
most fruitful areas for future research in accretion physics. Not only does an accurate 
treatment of the energetics allow for contact between theoretical models and observa- 
tions, but it also allows us to study the way in which the highly nonlinear temperature 
dependence of the resistivity, Hall and ambipolar diffusion terms influence the the non- 
linear behavior of the MRI. These results hold potentially significant promise for im- 
proving our understanding of the accretion process in protostellar accretion disks and 
outburst phenomena in dwarf nova disks. Even the modest question of how the inter- 
play between radiative cooling and turbulent heating influence the vertical structure of 
an accretion disk remains an open question. 
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